Projet à La Fin du Semestre
1 Modèle de régression linéaire
“Regression analysis is the hydrogen bomb of the statistics arsenal.”
Dans la première partie du projet, on va étudier comment à modéliser la rélation entre une variable aléatoire et plusieurs variables indépendantes. Cette méthode s’appelle la “régression linéaire.”
Pour le réaliser, d’abord, on va présenter un ensemble de données et l’analyser. Ensuite, on va trouver un modèle de régression linéaire convienant aux données. En fin, on va interpréter les résultats et tenter d’améliorer ce modèle.
1.1 Revue de littérature
L’analyse de régression peut être définie comme la recherche de la relation stochastique qui lie deux ou plusieurs variables. Son champ d’application recouvre de multiples domaines, parmi lequels on peut citer la physique, l’astronomie, la biologie, la chimie, la médecine, la géographie, la sociologie, l’histoire, l’économie, la lingustique et le droit. (Yadolah 2004, chap. 1)
La variable de résultat qu’on essaie d’expliquer s’appelle la variable expliquée, la variable dépendante, ou la réponse. Autres variables qu’on utilise à expliquer la réponse s’appellent les variables explicatives, les variables indépendantes, ou les prédicteurs.
1.2 Analyse exploratoire des données
On va construire un modèle de la régression linéaire pour le coût médical d’habitant dans quatre régions. Cet ensemble de données contient l’âge, le sexe, le BMI, le nombre d’enfants, les régions et une variable booléenne indiquée si un fume, qui sont indépendantes, et le coût médical étant unes variable dépendante.
Les données proviennent du forum de Kaggle (Kumar 2020), qui fournit un “terrain de jeu” pour des scientifiques de données. Dans l’intérêt de la simplicité, on ne utilise que l’âge, le sexe, le BMI et la variable de fumer comme ses variables indépendantes.
Regardons la sommaire des données générée par R :
insurance$age_group = cut(insurance$age, seq(10,70,10), c("10-20", "20-30", "30-40", "40-50", "50-60", ">60"), include.lowest=TRUE)
insurance$age_group = as.factor(insurance$age_group)
insurance$sex = as.factor(insurance$sex)
insurance$smoker = as.factor(insurance$smoker)
summary(insurance)## age sex bmi smoker charges
## Min. :18.00 female:662 Min. :15.96 no :1064 Min. : 1122
## 1st Qu.:27.00 male :676 1st Qu.:26.30 yes: 274 1st Qu.: 4740
## Median :39.00 Median :30.40 Median : 9382
## Mean :39.21 Mean :30.66 Mean :13270
## 3rd Qu.:51.00 3rd Qu.:34.69 3rd Qu.:16640
## Max. :64.00 Max. :53.13 Max. :63770
## age_group
## 10-20:166
## 20-30:278
## 30-40:257
## 40-50:281
## 50-60:265
## >60 : 91Une vérification rapide montre qu’il n’y a pas de valeur manquante dans son ensemble de données. Selon la sommaire, on voit que :
La moyenne d’âge des participant est plutôt haute : plus de 39 ans. On a divisé l’âge en groupe, créant une nouvelle variable explicative de “age_group.” Cependant, on va bénéficier cette variable seulement pour l’analyse des données; dans la partie de la modélisation plus tard, on va juste utiliser l’âge.
La ration entre les hommes et les femmes est approximative, mais pas celle du fumeur et du non-fumeur.
Les deux moyenne et médiane de l’indice de masse corporelle (BMI) sont plus grande que 30, que le WHO définit comme l’obésité. Cela peut être une signale d’un biais d’échantillonnage, bien qu’il soit trop tôt à conclure.
La réponse, “charges,” est bien rangée, de 1 122 à 63 770.
Maintenant, on va observer graphiquement quelques rélation entre les attributs dans ses données. Les codes sont adaptés de (Kassambara 2018).
Figure 1.1
library(ggplot2)
library(plotly)
cost_dist = ggplot(data = insurance, aes(charges)) +
geom_histogram(fill='steelblue',col='black', bins=20) +
geom_vline(xintercept = mean(insurance$charges), color = 'darkorange') +
geom_text(aes(x=mean(charges)+5000, label="\ncoût moyen", y=250), color="darkorange") +
labs(title = 'Fig 1.1. La Distribution du Coût Médical', y='Compte',x='Coût')
ggplotly(cost_dist + scale_color_brewer(palette="Dark2"))La distribution du coût médical est asymétrique à droite, avec la moyenne est environ 18270.
Figure 1.2
cost_by_smoker = ggplot(data = insurance, aes(charges, fill = smoker)) +
geom_histogram(bins=20, col='black') +
facet_wrap(~smoker) +
labs(title = 'Fig 1.2. Coût Médical par Fumeur', y='Compte',x='Coût')
ggplotly(cost_by_smoker + scale_color_brewer(palette="Dark2"))Généralement, les fumeurs doivent payer le coût médical beaucoup plus que les non-fumeurs. Un point notable est que dans son ensemble de données, la plupart de participants est non-fumeur (1 064/1 338). Si on trouve que la variable explicative de “smoker” est significative dans son modèle plus tard, on doit se rappeller ce point.
Figure 1.3
scatter = ggplot(insurance, aes(x=bmi, y=charges)) +
geom_point(color='steelblue') +
labs(title = 'Fig 1.3. Coût Médical par BMI', y='Coût', x='BMI')
ggplotly(scatter + scale_color_brewer(palette="Dark2"))Étonnamment, on ne peut pas voir une rélation claire entre le BMI et le coût médical des participants. Une commande rapide dans le logiciel de R le vérifie : la corrélation entre les deux variables est plutôt basse : 0.198341.
Figure 1.4
box = ggplot(data = insurance, aes(x=age_group, y=charges, fill=sex)) +
geom_boxplot() +
labs(title = 'Fig 1.4. Coût Médical par Sexe et Âge', y='Cost',x='Groupe d\'âge', color='Sexe')
box + scale_color_brewer(palette="Dark2")
Apparemment, le coût médical a une tendance à augmenter en fonction d’âge, mais le sexe ne diffère pas ce frais. La graphe au-dessus montre également qu’il y a beaucoup d’anomalie (outliers) dans son ensemble de données, qui fait la moyenne du coût à diffèrer bien de sa médiane (13 270 vs 9 382).
1.3 Modélisation
Après ayant recherché les données, maintenant, on essaie de les modéliser. Le modèle qu’on va inspecter est le modèle additif :
\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon \]
où :
\(X_1\): l’âge du participant.
\(X_2 = 1\) si le participant est masculin, sinon \(X_2 = 0\).
\(X_3\): l’indice de masse corporelle (BMI) du participant.
\(X_4 = 1\) si le participant fume, sinon \(X_4 = 0\).
Les hypothèses du modèle : le vecteur aléatoire \(\varepsilon\) suit une loi multinormale avec :
\[ E(\varepsilon) = 0 \\ Var(\varepsilon) = \sigma^2 \textbf{I}_n \]
On va faire appel au logiciel de R pour contruire son modèle :
model_add = lm(charges~age+sex+bmi+smoker, data = insurance)
model_add_summ = unlist(summary(model_add))
summary(model_add)##
## Call:
## lm(formula = charges ~ age + sex + bmi + smoker, data = insurance)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12364.7 -2972.2 -983.2 1475.8 29018.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11633.49 947.27 -12.281 <2e-16 ***
## age 259.45 11.94 21.727 <2e-16 ***
## sexmale -109.04 334.66 -0.326 0.745
## bmi 323.05 27.53 11.735 <2e-16 ***
## smokeryes 23833.87 414.19 57.544 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6094 on 1333 degrees of freedom
## Multiple R-squared: 0.7475, Adjusted R-squared: 0.7467
## F-statistic: 986.5 on 4 and 1333 DF, p-value: < 2.2e-16La droite de régression :
\[ y = -11633,49 + 259,45X_1 - 109,04X_2 + 323,05X_3 + 23833,87X_5 \]
C’est son modèle additif. Bien qu’il soit un peu simple et “sous-estimer” le problème, on va continuer à analyser ses résultats et essayer de l’améliorer en même temps.
1.4 Analyse des résultats
Signification de la régression
Le coefficient de détermination ajusté R2adj égale 0.7467396, qui indique une rélation linéaire plutôt forte entre la variable expliquée et les variables explicatives. On dit que ~75% de la variation observé des coût médical est expliquée par une rélation linéaire avec les variables explicatives.
Pour mieux comprendre la signification de la régression dans son modèle, on fait le test de Fisher avce l’hypothèse nulle suivant :
\[ \begin{aligned} & H_0: \beta_i = 0, \forall i \in \{1,2,3,4\} \\ & H_1: \exists i \in \{1,2,3,4\}: \beta_i \ne 0 \end{aligned} \]
On contruit le tableau d’ANOVA :
| Source de variation | Sommes des carrés | Degrés de liberté | Moyennes des carrés | Fobs |
|---|---|---|---|---|
| Régression | 146564944966 | 4 | 36641236241 | 986.5377 |
| Résiduelle | 49509276603 | 1333 | 37141243 | |
| Totale | 196074221568 | 1337 |
Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que la régression est significative, ou autrement dit, au moins une variable explicative a une forte rélation linéaire avec la réponse.
Tester sur chaque paramètre
Après on a testé la signification de la régression dans son modèle, on ensuite teste la signification de chaque variable explicative. On va faire 5 tests de Student, correspondant aux 5 variables différentes.
Rappeler son modèle est
\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon, \]
les tests sont présentés au-dessous.
Test 1
Les hypothèses :
\[ \begin{aligned} & H_0: \beta_0 = 0 \\ & H_1: \beta_0 \ne 0 \end{aligned} \]
Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle.
L’intervalle de confiance au niveau 95% pour \(\beta_0\) est :
conf_int[1,]## 2.5 % 97.5 %
## -13491.791 -9775.198Test 2
Les hypothèses :
\[ \begin{aligned} & H_0: \beta_1 = 0 \\ & H_1: \beta_1 \ne 0 \end{aligned} \]
Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_1\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.
L’intervalle de confiance au niveau 95% pour \(\beta_1\) est :
conf_int[2,]## 2.5 % 97.5 %
## 236.0266 282.8797Test 3
Les hypothèses :
\[ \begin{aligned} & H_0: \beta_2 = 0 \\ & H_1: \beta_2 \ne 0 \end{aligned} \]
Selon la sommaire du modèle, la valeur p-value = 0,745 > 0,05, donc on ne peut pas rejeter l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_2\) n’a pas d’effet dans son modèle. Ce résultat explique son observation dans la figure 1.4.
L’intervalle de confiance au niveau 95% pour \(\beta_2\) est :
conf_int[3,]## 2.5 % 97.5 %
## -765.5681 547.4858Test 4
Les hypothèses :
\[ \begin{aligned} & H_0: \beta_3 = 0 \\ & H_1: \beta_3 \ne 0 \end{aligned} \]
Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.
L’intervalle de confiance au niveau 95% pour \(\beta_3\) est :
conf_int[4,]## 2.5 % 97.5 %
## 269.0459 377.0563Test 5
Les hypothèses :
\[ \begin{aligned} & H_0: \beta_4 = 0 \\ & H_1: \beta_4 \ne 0 \end{aligned} \]
Selon la sommaire du modèle, la valeur p-value < 2.2e-16, donc on rejette l’hypothèse nulle au seuil \(\alpha\) = 5%. On conclut que \(\beta_0\) a un effet dans son modèle, ou autrement dit, l’âge a une rélation linéaire avec le coût médical.
L’intervalle de confiance au niveau 95% pour \(\beta_4\) est :
conf_int[5,]## 2.5 % 97.5 %
## 23021.34 24646.401.5 Diagnostic du modèle
Dans cette partie, on va diagnostiquer son modèle additif. C’est-à-dire, on va tester les hypothèses préséntés dans 1.2.
La régression linéaire fait plusieurs hypothèses sur les données, telles que :
Linéarité des données : la relation entre les prédicteurs et la réponse est censée être linéaire.
Normalité des résidus : les erreurs résiduelles sont supposées normalement distribuées.
Homogénéité de la variance des résidus : les résidus sont supposés avoir une variance constante (homoscédasticité).
Indépendance des résidus.
Les problèmes potentiels incluent :
Non-linéarité de la rélation entre réponse - prédicteurs.
Hétéroscédasticité : variance non-constante des résidus.
Présence de valeurs influentes dans les données (leverages et outliers).
Toutes ces hypothèses et problèmes potentiels peuvent être vérifiés en produisant des diagrammes de diagnostic visualisant les résidus. Le code pour le faire est adapté de (Rimal 2014).
Residual vs Fitted Plot
L’hypothèse de linéarité des données peut être vérifiée en inspectant la figure de Residual vs Fitted.
Idéalement, la ligne orange doit être proche de la droite rouge; dans ce cas la rélation entre \(Y\) et \(X\) est linéaire. Mais la graphe au-dessus indique une rélation non-linéaire dans les données.
Normal Q-Q
L’hypothèse de normalité des résidus peut être vérifiée en inspectant la figure de Normal Q-Q.
Si les résidus suivent une loi normale, ils ont être approximatifs la droite orange. Dans son modèle, il semble qu’il est en opposé. On peut le tester en faisant le test de Shapiro-Wilk.
\[ \begin{aligned} & H_0: \varepsilon \text{ suit une loi normale} \\ & H_1: \varepsilon \text{ ne suit pas une loi normale} \end{aligned} \]
Le logiciel de R a déjà une commande pour le faire :
res = resid(model_add)
shapiro.test(res)##
## Shapiro-Wilk normality test
##
## data: res
## W = 0.90294, p-value < 2.2e-16Avec le valuer p-value < 2.2e-16, on rejette l’hypothèse nulle au seil \(\alpha\) = 5%; il s’agit du non-normalité des résidus.
Scale-Location
L’hypothèse d’homogénéité de la variance des résidus peut être vérifiée en inspectant la figure de Scale-Location.
C’est bien si on voit une ligne horizontale avec des points également répartis. Dans son exemple, ce n’est pas le cas. La variance a tendance à augmenter avec les valeurs “fitted” de la réponse, qui insinue l’hétéroscédasticité. On peut le tester en faisant le test de Breusch-Pagan.
\[ \begin{aligned} & H_0: \text{Les résidus sont distribués avec une variance constante} \\ & H_1: \text{Les résidus ne sont pas distribués avec une variance constante} \end{aligned} \]
Le logiciel de R a déjà une commande pour le faire :
library(lmtest)
bptest(model_add)##
## studentized Breusch-Pagan test
##
## data: model_add
## BP = 113.48, df = 4, p-value < 2.2e-16Avec le valuer p-value < 2.2e-16, on rejette l’hypothèse nulle au seil \(\alpha\) = 5%; cela indique l’hétéroscédasticité des résidus.
Residual vs Leverage
On peut trouver des outliers dans les données en inspectant la figure de Scale-Location.
On voit qu’il y a beacoup des points de donnée qui dépasse 3 écart-types de la droite de régression. Cela signifie que son ensemble de données se composent des outliers, qui affecte négativement la performance de son modèle.
1.6 Améliorations possibles
Transformation de données
La plupart des problèmes qu’on a rencontré dans la section précédente peuvent être résolus en appliquant une transformations (logarithme, racine carrée…) à la variable de résultat \(Y\). On va faire la transformation de Box-Cox avec \(\lambda=0,2\) dans ce cas :
model_bc = lm(
((charges)^0.2 - 1) / 0.2 ~ age + sex + bmi + smoker,
data=insurance)
summary(model_bc)##
## Call:
## lm(formula = ((charges)^0.2 - 1)/0.2 ~ age + sex + bmi + smoker,
## data = insurance)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2792 -1.4073 -0.3976 0.5023 13.0296
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.937013 0.434849 32.050 < 2e-16 ***
## age 0.200296 0.005482 36.538 < 2e-16 ***
## sexmale -0.342500 0.153630 -2.229 0.026 *
## bmi 0.086518 0.012637 6.846 1.15e-11 ***
## smokeryes 10.269389 0.190135 54.011 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.798 on 1333 degrees of freedom
## Multiple R-squared: 0.7627, Adjusted R-squared: 0.7619
## F-statistic: 1071 on 4 and 1333 DF, p-value: < 2.2e-16Bien que le R2adj soit mieux que le premier modèle, les tests de Shapiro-Wilk et de Breusch-Pagan donnent les conclusions similaires : il n’y a pas de normalité et de homogénéité des résidus.
shapiro.test(resid(model_bc))##
## Shapiro-Wilk normality test
##
## data: resid(model_bc)
## W = 0.85154, p-value < 2.2e-16bptest(model_bc)##
## studentized Breusch-Pagan test
##
## data: model_bc
## BP = 53.586, df = 4, p-value = 6.425e-11Note : des détails de la méthodologie de Box-Cox sont hors sujet du cours, faire référence à (Dalpiaz 2018, chap. 14) pour la lecture supplémentaire.
Conclusion : la transformation de la variable de résultat n’est pas d’effet significatif.
Interaction entre les prédicteurs - Multicolinéarité
Dans le commencement de ce rapport, on a mentionné que le modèle additif est probablement une sous-estimation du problème. En fait, c’est exact. Dans son modèle additif, le coût médical peut être différent en moyenne entre smoker et non-smoker pour le même âge; mais le changement du coût pour une augmentation d’âge est égale pour les deux groupes.
Cela vient du fait que les interactions entre les prédicteurs ne sont pas inclus dans son modèle. Cettes interactions sont très importantes, surtout si on a des variables qualitatives comme smoker et sex.
Inspester le modèle avec l’interaction est plutôt complexe, sans parler du problème de multicolinéarité. Car, une fois encore, dans l’intérêt de la simplicité, on ne va pas le présenter ici et le laisser pour la lecture supplémentaire.
1.7 Conclusion
Après beaucoup d’effort de modéliser les données, son modèle est :
\[ Y'(X) = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_4X_4 + \varepsilon \]
où :
\(X_1\): l’âge du participant.
\(X_2 = 1\) si le participant est masculin, sinon \(X_2 = 0\).
\(X_3\): l’indice de masse corporelle (BMI) du participant.
\(X_4 = 1\) si le participant fume, sinon \(X_4 = 0\).
\(Y' = log(Y)\), où \(Y\) est le coût médical.
La droite de régression :
\[ Y'(X) = 13,937 + 0,2X_1 - 0,343X_2 + 0.087X_3 + 10.269X_4 \]
Les résultats peuvent être décevant, mais quelquefois, on doit accepter que sa technique ne peut pas résoudre le problème sous la main, surtout en réalité. C’est pourquoi on a toujours besoin d’apprendre de nouvelles choses et de les ajouter à sa boîte à outils !
2 Choix du modèle
“All models are wrong, but some are useful.”
Lorsqu’on fait l’analyse de régression, on se rend compte souvent qu’on doit chercher dans un trop grand espace de variables explicatives. C’est difficile, parfois impossible, à rechercher exhaustive, c’est-à-dire à essayer toutes les combinaisons des prédicteurs. C’est sans parler des écueils potientiels présentés dans la section précédente. Car, dans cette partie, on va démontrer deux meilleurs procédure pour établir un modèle de régression, à savoir la méthode de stepwise et la méthode de stagewise.
Note : Critères de choix
On a utilisé un critère de choix dans la première section : c’est la valeur du R2 ajusté. En fait, le critère du R2 ajusté et du R2 sont les plus utilisés, bien qu’il y aie cependant des critiques. Il existe également plusieurs critères pour selectionner le meilleur modèle :
Critère du Cp de Mallows : l’idée de ce critère est de choisir un modèle pour lequel la somme de erreurs quadratiques moyenne est minimale. On choisira la valeur du coefficient Cp la plus proche de nombre de paramètres \(p\) dans le modèle.
Critère d’information d’Akaike (AIC) : se compose la vraisemblance (the likelihood) qui mesure “la qualité de l’ajustement,” et la pénalité (the penalty) qui est fonction de la taille du modèle. On choisira la valeur d’AIC la plus petite.
Critère d’information bayésien (BIC) : quantifie le compromis entre un modèle qui s’ajuste bien et le nombre de paramètres du modèle, mais pour une taille d’échantillon raisonnable, il choisit généralement un modèle plus petit que AIC.
2.1 Exercise 1 : Mensurations
On s’intéresse au lien éventuel entre le poids d’un homme et divers caractéristiques physiques. Son ensemble de données se compose 22 hommes en bonne santé âgés de 16 à 30 ans.
On va utiliser la procédure de stepwise avec le critère d’AIC pour trouver un modèle de régression pour ces données.
Procédure de stepwise
La procédure stepwise propose après l’introduction d’une nouvelle variable dans le modèle :
Réexaminer les tests de Student pour chaque variable explicative anciennement admise dans le modèlel,
Après réexamen, si des variables ne sont plus significatives, alors retirer du modèle la moins significative d’entre elles.
Démonstration
On peut faire appel au logiciel de R à effectuer la procédure de stepwise :
model1_start = lm(Y~1, data = measurement)
model1_both_aic = step(
model1_start,
scope = Y ~ X1+X2+X3+X4+X5+X6+X7+X8+X9+X10,
direction = 'both'
)## Start: AIC=106.34
## Y ~ 1
##
## Df Sum of Sq RSS AIC
## + X6 1 2113.64 410.50 68.380
## + X1 1 2038.88 485.27 72.060
## + X5 1 1852.89 671.26 79.198
## + X9 1 1790.72 733.43 81.147
## + X8 1 1747.28 776.87 82.413
## + X4 1 1648.10 876.05 85.056
## + X3 1 1519.92 1004.23 88.061
## + X2 1 1333.64 1190.50 91.804
## + X7 1 610.86 1913.29 102.242
## <none> 2524.15 106.338
## + X10 1 153.59 2370.56 106.956
##
## Step: AIC=68.38
## Y ~ X6
##
## Df Sum of Sq RSS AIC
## + X1 1 249.97 160.54 49.724
## + X8 1 194.80 215.71 56.223
## + X5 1 172.36 238.14 58.400
## + X9 1 151.74 258.76 60.227
## + X7 1 83.70 326.80 65.363
## + X4 1 82.56 327.95 65.440
## + X3 1 76.14 334.37 65.866
## + X2 1 63.22 347.29 66.700
## <none> 410.50 68.380
## + X10 1 2.11 408.39 70.266
## - X6 1 2113.64 2524.15 106.338
##
## Step: AIC=49.72
## Y ~ X6 + X1
##
## Df Sum of Sq RSS AIC
## + X7 1 48.39 112.15 43.833
## + X8 1 28.85 131.68 47.366
## + X9 1 23.45 137.08 48.250
## <none> 160.54 49.724
## + X10 1 11.59 148.94 50.076
## + X5 1 8.94 151.60 50.464
## + X2 1 7.30 153.24 50.701
## + X3 1 1.13 159.41 51.569
## + X4 1 0.03 160.50 51.720
## - X1 1 249.97 410.50 68.380
## - X6 1 324.73 485.27 72.060
##
## Step: AIC=43.83
## Y ~ X6 + X1 + X7
##
## Df Sum of Sq RSS AIC
## + X9 1 26.174 85.97 39.986
## + X8 1 17.481 94.67 42.105
## + X10 1 13.685 98.46 42.970
## <none> 112.15 43.833
## + X3 1 5.196 106.95 44.790
## + X4 1 0.803 111.35 45.675
## + X5 1 0.272 111.88 45.780
## + X2 1 0.120 112.03 45.810
## - X7 1 48.387 160.54 49.724
## - X1 1 214.655 326.80 65.363
## - X6 1 285.702 397.85 69.691
##
## Step: AIC=39.99
## Y ~ X6 + X1 + X7 + X9
##
## Df Sum of Sq RSS AIC
## + X10 1 12.019 73.955 38.673
## + X8 1 8.480 77.494 39.701
## <none> 85.974 39.986
## + X4 1 1.537 84.437 41.589
## + X5 1 1.119 84.855 41.698
## + X3 1 0.160 85.814 41.945
## + X2 1 0.061 85.913 41.970
## - X9 1 26.174 112.149 43.833
## - X7 1 51.107 137.081 48.250
## - X1 1 97.997 183.971 54.722
## - X6 1 194.495 280.469 63.999
##
## Step: AIC=38.67
## Y ~ X6 + X1 + X7 + X9 + X10
##
## Df Sum of Sq RSS AIC
## + X8 1 9.233 64.722 37.739
## <none> 73.955 38.673
## + X3 1 3.374 70.581 39.646
## - X10 1 12.019 85.974 39.986
## + X2 1 1.191 72.764 40.316
## + X4 1 0.794 73.162 40.436
## + X5 1 0.319 73.636 40.578
## - X9 1 24.508 98.464 42.970
## - X7 1 53.031 126.986 48.567
## - X1 1 105.028 178.983 56.117
## - X6 1 203.165 277.120 65.735
##
## Step: AIC=37.74
## Y ~ X6 + X1 + X7 + X9 + X10 + X8
##
## Df Sum of Sq RSS AIC
## <none> 64.722 37.739
## + X3 1 4.748 59.974 38.063
## - X8 1 9.233 73.955 38.673
## + X2 1 1.593 63.129 39.191
## + X5 1 1.350 63.372 39.276
## - X10 1 12.772 77.494 39.701
## + X4 1 0.004 64.718 39.738
## - X9 1 15.559 80.281 40.479
## - X7 1 42.815 107.538 46.910
## - X1 1 59.005 123.727 49.995
## - X6 1 196.988 261.710 66.476Le résultat de la procédure est le modèle :
\[ Y = \beta_0 +\beta_1X_6 + \beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_8 + \varepsilon \]
L’explication étape par étape :
Étape 1 : on commence avec le modèle sans variable explicative. On essaie d’ajouter une variable parmi les dix variables à son modèle. Comme le modèle linéaire avec la variable \(X_6\) a la meilleur statistique d’AIC , on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 + \varepsilon\).
Étape 2 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_1\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \varepsilon\).
Étape 3 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_7\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \varepsilon\).
Étape 4 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_9\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \varepsilon\).
Étape 5 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_10\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \varepsilon\).
Étape 6 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. L’addtion de variable \(X_{8}\) a la meilleur statistique d’AIC , donc on le sélectionne et continue.
Le modèle obtenu après cette étape : \(Y = \beta_0 +\beta_1X_6 +\beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_{8} + \varepsilon\).
Étape 7 : on essaie d’ajouter une nouvelle variable à son modèle en faisant le test de Student pour chaque variable anciennement admise dans le modèle. Maintenant, aucune de modification du modèle précédant peut amélioré sa performance, donc on arrête et conclut.
Le modèle final de son choix :
\[ Y = \beta_0 +\beta_1X_6 + \beta_2X_1 + \beta_3X_7 + \beta_4X_9 + \beta_5X_{10} + \beta_6X_8 + \varepsilon \]
La droite de régression :
\[ Y = -79,73 + 0,66X_6 + 1,79X_1 + 0,25X_7 + 0,43X_9 - 0,66X_{10} + 0,51X_8 + \varepsilon \]
Conclusion
On a établi le meilleur modèle pour son ensemble de données par la procédure de stepwise. Alors que cette procédure ne sélectionne pas nécessairement le meilleur modèle absolu, elle reste généralement un modèle acceptable.
2.2 Exercise 2 : Taux d’accidents
En inspectant 39 observations faites sur des troncons d’autoroute, on veut trouver un modèle de régression linéaire pour expliquer le taux d’accidents dans l’état du Minnesota.
On va utiliser la procédure de stagewise pour trouver un modèle de régression pour ces données.
Procédure de stagewise
Cette méthode se déroule de la façon suivante :
Effectuer la régression avec la variable la plus corrélée avec Y.
Calculer les résidus obtenus avec cette régression.
Considérer ensuite ces résidus comme une nouvelle variable dépendante que l’on veut expliquer à l’aide des variables explicatives restantes.
Démonstration
Dans cet exercise, on considère une corrélation plus petite de 0,2 est insignifiante.
- Étape 1 : on recherche la variable explicative \(X_i\) la plus corrélée avec \(Y\). On a
## Cor.vs.Y
## x_i,9 0.752025471
## x_i,4 -0.680983625
## x_i,8 0.564479507
## x_i,3 -0.512522209
## x_i,1 -0.428100278
## x_i,6 -0.386907190
## x_i,13 0.337848301
## x_i,11 -0.207610359
## x_i,12 -0.161539823
## x_i,7 -0.160385294
## x_i,10 -0.032978589
## x_i,2 -0.028569805
## x_i,5 -0.005619311Parce que \(X_9\) est la plus corrélée avec \(Y\), on fait la régression linéaire avec lui comme un seule variable explicative. On obtient la valeur estimée \(\hat y_i\) à partir de cette esquation :
\[ \hat y_i = 1,98 + 0,16x_{i, 9} \]
Étape 2 : on calcule les résidus \(e_{i,1} = y_i - \hat y_i\) du modèle précédent.
Étape 3 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,1}\) au-dessus.
## Cor.vs.Residus
## x_i,1 -0.41525386
## x_i,3 -0.36653065
## x_i,8 0.28742153
## x_i,4 -0.25558207
## x_i,2 0.21220066
## x_i,10 0.18816837
## x_i,7 0.11808444
## x_i,6 -0.10215374
## x_i,11 0.07597809
## x_i,13 -0.07328953
## x_i,5 0.03940902
## x_i,12 0.01466631Parce que \(X_1\) est la plus corrélée avec \(e_{i,1}\) calculé à l’étape 2, on fait la régression linéaire entre ces résidus et \(X_1\) comme un seule variable explicative. On obtient la valeur estimée \(e_{i,1}\) à partir de cette esquation :
\[ \hat e_{i,1} = 0,88 - 0,07x_{i, 1} \]
Étape 4 : on calcule les résidus \(e_{i,2} = e_{i,1} - \hat e_{i,1}\) du modèle précédent.
Étape 5 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,2}\) au-dessus.
## Cor.vs.Residus
## x_i,4 -0.22180561
## x_i,3 -0.20396676
## x_i,8 0.18163358
## x_i,6 -0.17319831
## x_i,5 -0.09634904
## x_i,10 0.08760933
## x_i,2 0.08175689
## x_i,12 -0.03793432
## x_i,11 0.03468900
## x_i,7 0.03091557
## x_i,13 -0.01328866Parce que \(X_4\) est la plus corrélée avec \(e_{i,1}\) calculé à l’étape 4, on fait la régression linéaire entre ces résidus et \(X_4\) comme un seule variable explicative. On obtient la valeur estimée \(e_{i,2}\) à partir de cette esquation :
\[ \hat e_{i,2} = 2,48 -0,05x_{i, 4} \]
Étape 6 : on calcule les résidus \(e_{i,3} = e_{i,2} - \hat e_{i,2}\) du modèle précédent.
Étape 7 : on recherche la variable explicative \(X_i\) la plus corrélée parmi les variables explicatives restantes avec \(e_{i,3}\) au-dessus.
## Cor.vs.Residus
## x_i,10 0.15001818
## x_i,3 -0.14180361
## x_i,11 0.14133763
## x_i,2 0.13938413
## x_i,13 -0.11006480
## x_i,7 0.10945059
## x_i,8 0.09296005
## x_i,5 -0.07636050
## x_i,12 -0.02879767
## x_i,6 -0.02089275On se rend qu’il n’y a pas de variable explicative étant corrélée significativement avec \(e_{i,3}\). La procédure de stagewise s’arrête donc là.
Son équation final est l’addition de deux équations obtenues aux étapes 1, 3 et 5 :
\[ \begin{eqnarray} y_i &=& \hat y_i + e_{i,1} \\ &=& 1,98 + 0,16x_{i, 9} + \hat e_{i,1} + e_{i,2} \\ &=& 1,98 + 0,16x_{i, 9} + 0,88 - 0,07x_{i, 1} + e_{i,2} \\ &=& 2,86 + 0,16x_{i, 9} - 0,07x_{i, 1} + \hat e_{i,2} + e_{i,3} \\ &=& 2,86 + 0,16x_{i, 9} - 0,07x_{i, 1} + 2,48 -0,05x_{i, 4} + e_{i,3} \\ &=& 5,34 + 0,16x_{i, 9} - 0,07x_{i, 1} - 0,05x_{i, 4} + e_{i,3} \\ \end{eqnarray} \]
Le modèle équivalent est :
\[ Y = 5,34 + 0,16X_{9} - 0,07X_{1} - 0,05X_{4} + \varepsilon \]
Conclusion
Grâce à la procédure de stagewise, on peut établir un modèle de régression linéaire pour son ensemble de données. Cependant l’estimation par les moindres carrés fournit une prédiction globale meilleure que la régression stagewise, elle offre de bons résultats quand on suspecte qu’il y aie un problème de multicolinéarité.
3 Référence
Réalisé par Groupe 1
Khang Pham | Phat Ngo | Chau Ho